Work distributions in the T = Random Field Ising Model 
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We perform a numerical study of the three-dimensional Random Field Ising Model at T = 0. We 
compare work distributions along metastable trajectories obtained with the single-spin flip dynamics 
with the distribution of the internal energy change along equilibrium trajectories. The goal is to 
investigate the possibility of extending the Crooks fluctuation theorem [l[ to zero temperature when, 
instead of the standard ensemble statistics, one considers the ensemble generated by the quenched 
disorder. We show that a simple extension of Crooks fails close to the disordered induced equilibrium 
phase transition due to the fact that work and internal energy distributions are very asymmetric. 



PACS numbers: 75.60.Ej, 75.50.Lk, 81.30.Kf, 81.40.Jj 



£ ' I. INTRODUCTION 

C ' 

C/3 . The T = Random Field Ising Model is a prototype model for the study of collective phenomena in disordered 
systems. Although it neglects thermal fluctuations, it contains essential competitions between the quenched disorder, 
the ferromagnetic interaction and the external applied field. The model can be numerically studied from two different 
points of view: on the one hand, the exact ground state calculation 0, SIHUI provides an approach to the equilibrium 
phase diagram. On the other the use of a local relaxation dynamics based on single spin-flips provides a good framework 
for the understanding of avalanche dynamics and hysteresis 0, 0], which is closer to experimental observations. In 
this sense, the model is a good workbench for the comparison of equilibrium and out-of-equilibrium trajectories. 

A number of non-equilibrium work theorems @ have received a lot of attention in the last 10 years, particularly 
after the work of Jarzynski in 1997 HEH. These theorems relate in different ways the distribution of work performed 
on a system which is driven out of equilibrium to some equilibrium thermodynamic properties. One example is the 
original Jarzinsky's equality: (e _/3M/ ) = e~ l3AF , where W is the work performed on the system that has been driven 
(out-of-equilibrium) by varying an external control parameter H changing from H(0) to H(l), and AF is the free 
energy difference between two states and 1 that correspond to the equilibrium states at H(0) and H (1). The average 
(•) should be understood as obtained after many repetitions of the driving process. The system is assumed to be in 
, contact with a heat bath at temperature which generates an statistical ensemble of copies of the system which is 
t^J" ' the source of work fluctuations. Most of these theorems have been experimentally verified. [Til. [l2| 

The goal of this paper is to investigate how such kind of theorems can be extended to systems at T = 0. In 
t—( ■ such a case, although equilibrium states and out-of-equilibrium trajectories are well defined, there are no thermal 
fluctuations. A priori it seems that there is no statistical ensemble over which one can define probability distributions 
or averages. The idea we want to test is whether or not, for systems with quenched disorder, the thermal ensemble 
can be substituted by the ensemble of different realizations of disorder and still some work theorems can be applied. 
• i-H . In order to perform the T — ► limit we choose a different but related work theorem which is Crooks fluctuations 
theorem [l|, [l3| . The advantage will be that it allows to derive an equality that can be extrapolated to the T — > 



^ , limit. This theorem can be written as: 
C3 
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In this case Pp(W) and Pr(—W) are the probability densities corresponding to the out-of-equilibrium work performed 
on a system driven forward from H(0) to -ff(l) and reversely from H(l) to H(0). The Crooks fluctuation theorem 
has been formulated under several assumptions: the driven systems must be finite, classical and coupled to a thermal 
bath. Moreover, the dynamics should be stochastic, Markovian and reversible and the entropy production should be 
odd under time reversal. Some of these assumptions are clearly not accomplished at the work at hand. For instance 
the thermal bath is at T — and the dynamics is deterministic. Nevertheless we have been a little speculative and 
investigated the possibility of extending the Crooks fluctuation theorem to T — for systems with quenched disorder. 

From Eq. one can derive that the value W* for which the two probability densities are equal, i.e. Pf(W*) = 
Pr{— W*) satisfies W* = AF. This result is particularly suitable to be investigated at T — > 0. Note that in the 
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low temperature limit, the only possibility to have a crossing point of the probability densities at W* is that W* = 
AF(T — *• 0) = AU, so that the diverging behaviour of (3 might be cancelled to get a finite limit P F (W*) / P R {-W*) — » 
1. 



II. MODEL 

For our investigation we consider the Random Field Ising Model. A set of N spin variables Si = ±1 with i = 1, • • • N 
is defined on a regular 3D cubic lattice with linear size L (N = Lx Lx L). We are interested in following the response 
of the order parameter (magnetization) Al = $^ i= i as a function of the external driving field H . The Hamiltonian 
(magnetic enthalpy) of the model is defined as: 

H = U - HM (2) 

N 

U = -J2 S i S i - hiSi ( 3 ) 

ij i=l 

where the first sum in Eq.[3] extending over nearest neigbours pairs, accounts for a ferromagnetic exchange interaction, 
the second term in Eq. [3] includes the interaction with the quenched disorder and the second term in Eq. [2] accounts 
for the interaction with the external field H that will be used as the driving parameter. The local fields hi are 
independent random variables, quenched, Gaussian distributed with zero mean and standard deviation a. All along 
the paper we will use small letters (u = j$,m = 4£ . . .) to refer to intensive magnitudes. 

The equilibrium properties of this model can be studied at T = 0. For a given realization of the random fields 
{hi} and a given value of H, the exact ground state can be found in a polynomial time [14j . Moreover optimized 
algorithms can also be used to obtain the full sequence of ground-sates as a function of H between the two saturated 
states {Si = —1} at H = — oo and {Si — 1} at H = oo [la. In the thermodynamic limit the 3D model exhibits a first 
order phase transition at H = and for a < cr^ q = 2.27 @, H|. Although finite systems do not exhibit a true phase 
transition, there is a region of a and H where the correlation length is of the same order as the system size. This 
leads to collective effects and correlations involving all the spins of the system. 

Concerning the non-equilibrium trajectories, the model has been intensively studied by using the so called single- 
spin flip metastable dynamics at T — which can be understood as a T — > Glauber dynamics. It consists in 
adiabatically sweeping the external field and relaxing individual spins according to their local energy, i.e. spins should 
align with the local field 

Fi = H + hi + S k (4) 

k 

where the sum extends over the neigbours of spin i. Note that the reversal of a spin at a given field H may induce 
the reversal of some of its neigbours at the same field H . Such collective events are the so called avalanches that end 
when all the spins are stable. Only after the avalanches have finished, the external field is varied again. It has been 
shown that, during a monotonous driving of the external field, no reverse spin-flips may occur and, moreover, the 
model exhibits the so called abelian property: i.e. the final states do not depend on the order in which the unstable 
spins are flipped. In the thermodynamic limit, the system with this metastable dynamics also shows a critical point 
at cr™ et = 2.21 @, [IB], below which the metastable trajector y ex hibits a discontinuity. 

Note that avalanches are the source of energy dissipation [l7l . Il8j . If we think about an increasing field trajectory, 
the triggering spin of each avalanche flips at the field H that corresponds to Fi = 0. This spin, therefore, flips without 
energy loss AH = AU — HAM = 0. But the subsequent unstable spins flip (at a fixed external field H) with an 
external force Fi > giving raise to an energy loss (associated to each individual spin flip) Q = —2Fi < with 
Q = AU — W. For the decreasing field branches, Fi < but the energy loss is then Q = 2Fi < since the spins flip 
from 1 to —1. 

Note that in this discussion we are using the standard definition of work that is not a state function. It is computed 
as a sum over the (k = 1, 2 . . .) spins that flip along the trajectory 

W = J2 H k^M k , (5) 

k 

where AM^ — 2. Recently there has been a discussion 0, [53, Hi], H2| on which is the most suitable definition of 
work to be used in such non-equilibrium work theorems. Already from the initial Jarzynski works [9(, it was proposed 
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that if the system is driven by controlling the external field, the convenient definition of work is the integral over 
the trajectory V — J MdH. Note that along a metastable trajectory , the two definitions of work are related by 
W — A{HM) — V. Without going into the discussion, we will numerically test here the two possibilities. In Sec. IIIII 
we will test whether the crossing point W* of the histograms Pp(W) and Pr(—W) satisfies 

W* = (AU) (6) 

and, in Sec. IIVI we will test whether the crossing point V* of the histograms Pf(V) and Pr(—V) satisfies the 
corresponding equation: 

V* = (AU - A(HM)) (7) 

The test of the two hypothesis will require slightly different strategies (as indicated in Fig. [J) since V can not be 
computed for trajectories starting at saturation (H —> ±oc). In the first strategy (Sec. IHI|) . we will compute forward 




FIG. 1: Schematic representation of the two different strategies that have been used, (a) corresponds to the first strategy in 
which the starting point is the saturated state (m = —1,H = — oo), and (b) corresponds to the second strategy, starting from 
the equilibrium state at a finite field Ho- Dashed lines correspond to metastable trajectories and the continuum line indicates 
the equilibrium states. 



trajectories starting from the state 0, saturated with m = —I at H = — oo and sweep the field adiabatically from — oo 
to H\ until the metastable state A is reached. Then we will compute the equilibrium state 1 corresponding to the field 
Hi and perform the reverse trajectory from 1 to with the metastable dynamics. We will compute Ait = (Ui — Uq) /N 
and the following works (using the 'standard' definition): 

w i if 1 i i 

w ^x = = ~r J o HdM +n] a HdM =nJ HdM + N H{Ml - Ma) (8) 



N N 
Wi^o 1 







Wl ^° = = NJi HdM (9) 

where Ma and Mi are the magnetizations of states A and 1. The integrals are computed along the trajectories 
schematically represented in Fig. [Ua) usm g Eq. [5l Note also that the second integral in Eq. [5] can be computed since 
the process A — * 1 takes place at constant field Hi. 

In the second strategy (Sec. IIV|) . we will start from a computed equilibrium state at Hq, perform a metastable 
trajectory until the state A is reached at H\. Then we will compute the equilibrium state 1 at Hi and perform a 
metastable trajectory driving back until reaching the state B at Hq. In this case we will use the alternative definition 
of work: 

t/„ . TA - 1 r A 

MdH, (10) 

B 

MdH. (11) 
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The integrals are computed along the trajectories schematically represented in Fig.QJb) using V = ^2 k M^AH^ where 
k accounts for the sequence of interavalanche field increments AH^ occuring at constant M^ along the trajectory. 
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FIG. 2: (Color on line) Example of histograms corresponding to Pf(wo^i), Pfl(-iui^o) and P(Au). The arrow pointing 
downwards indicates the position of the maximum of P(Au) and the arrow pointing upwards the average value (Au). The 
example corresponds to a system with a = 4, Hi — —0.35 and L — 12. Histograms have been obtained by cumulating data 
corresponding to 4 x 10 5 realizations of disorder. 



Fig. 2 shows histograms corresponding to the distributions of wq^i, Au and — wi^q for er = 4, Hi = —0.35 
and a system size L = 12. The computed histograms give an estimation of the corresponding probability densities 
Pf(wq^i), Pr{— w\^q) and P(Au). Four important points should be realized from this example: (i) For this value of 
a and field Hi, the densities look symmetric and Gaussian. This will not be the case when Hi and a are close to the 
disorder induced phase transition at a^ q and Hi — 0. (ii) Second, the distribution of Au has a width similar to those 
of the out-of-equilibrium works. This is different from what happens in the analysis of Crooks fluctuation theorem 
at finite T. Typically, equilibrium thermal fluctuations are much smaller than work fluctuations, (iii) For increasing 
systems size the histograms become narrower and then the crossing points are harder to locate (iv) Finally, note that 
for this case the hypothesis that we are testing is fulfilled: the peak in P(Au) as well as the average (Au) coincide with 
the crossing point w* of Pr(wq^i) and Pr{— ioi_>o)- Fig.[3]shows histograms corresponding to P{wq— >i), P(Au) and 



FIG. 3: (Color on line) Example of histograms corresponding to Pf(wo^i), Pr(-wi^o) and P(Au). The arrow pointing 
downwards indicates the position of the maximum of P(Au) and the arrow pointing upwards the average value (Au). The 
example corresponds to a system with a = 3, Hi =0.1 and L = 12. The histograms have been computed by cumulating data 
corresponding to 1.6 x 10 6 realizations of disorder. The inset shows a detailed view of the crossing point of the Pf(uio-.i) and 
Pr(— wi->o) histograms. 

P(— lui-fo) f° r a — 3, Hi = 0.1 and a system size L — 12. Note than in this case, the work distributions as well as 
the energy distribution are very asymmetric. This is because although in this case we expect Mi > since Hi > 0, a 
certain non vanishing fraction of the equilibrium states still has Mi < 0, as schematically indicated in FigUJ In other 



III. FIRST STRATEGY 
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FIG. 4: Schematic representation of the case a < a^ q and H\ > 0. Dashed lines correspond to metastable trajectories, contin- 
uous lines to the equilibrium trajectory and the dashed-dotted line to the non typical equilibrium states that are responsible 
for the high asymmetry of the distribution Pr(— wi-to)- 

words, for certain realizations of disorder, the state 1 turns out to be non-typical and previous to the equilibrium 
transition from M < to M > 0. Consequently, the distribution of reverse works uq^o widens enormously. As can 
be seen the proposed equality w* = Au clearly fails in this case. The relative error is \w* — (Au)|/|(Au)| ~ 0.25. 

We have performed a detailed study of the deviation of w* from (Au) for different values of a, Hi and L. Examples of 
the results are presented in Fig. [5] Note that the data corresponding to the crossing point w* exhibits increasing error 
bars with increasing Hi. This is because the histograms of Pf(wo^i), Pr(— wi->o) become more and more separated 
and the finding of the intersection requires more and more statistics. (This problem is accentuated for larger values 
of Hi and for larger system sizes). One can conclude therefore than the proposed extrapolation of Crooks theorem, 
given by the equality w* = (Au) fails in the region of a and H\ where the system exhibits a collective behaviour 
(i.e. the correlation length is similar to the system size) because of the proximity of the equilibrium phase transition. 
Consistently, we observe in Fig. Efb) (comparing the data corresponding to L — 12 and L — 18 at Hi ~ 0) that the 
region of breakdown becomes smaller when the system size L is larger. As expected, increasing the system size with 
fixed correlation length decreases the collective effects. 

IV. SECOND STRATEGY 

Fig. [5]shows histograms corresponding to the distributions of vq—,i and — i>i^o for the case L — 12, a = 3, Hi = — 10 
and H2 = —0.4. As can be seen the intersection point v* agrees very well with the average value of Au — A(Hm). In 
this case we have also studied if this agreement is equally good in other regions of the phase diagram. Fig [7] shows 
a comparison of (Au) and v* + (A(Hm)}. On the right hand side, the value of v* has not been obtained by a direct 
location of the intersection of the histograms but locating the mid point between the positions of the maxima of the 
histograms corresponding to Pf(uo_>i) and Pr{— t?i— >o) which are far separated but very symmetric. As can be seen, 
as occurs with strategy 1, the agreement v* ~ (Au — A(Hm)) fails when Hi approaches the phase transition region. 

V. SUMMARY AND CONCLUSION 

We have investigated the possibility of extending some work fluctuation theorems to T = for systems with 
quenched disorder. In this case, thermal averages should be substituted by disorder averages. We have proposed 
an hypothesis based on the extrapolation of Crooks theorem to the T — > limit and we have numerically tested its 
validity. The hypothesis is that if one considers the distributions of non-equilibrium work corresponding to a forward 
and backwards trajectory, the crossing point w* of Pf(iuo-»i) and Pr(—wi->o), is equal to the average of the internal 
energy difference Au = u\ — u%. The investigation has been done by considering two strategies, based on the two 
possible definitions of work that have been discussed in the literature: w — J Hdm and v = J mdH. 

The reported numerical investigation indicates that, for both strategies, the formulated hypothesis is valid when 
the system does not behave collectively. Thus, far from the equilibrium critical point, where correlations are small 
compared to system size, the distribution p(Au) is very symmetric and the average (equilibrium) value (Au) can be 
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FIG. 5: (Color on line) Comparison of the crossing point of the histograms w* (□) with (Am) (O) as a function of Hi, 
corresponding to (a) a — 4 and (b) a = 3 for different values of L as indicated by the legend. For the case L = 12 and a — 3 
the dashed line indicated the position of the maximum of P(Au). Data for L — 18 (L = 24) has been shifted 0.3 (0.6) units 
along the vertical axis in order to clarify the picture. 

obtained from the crossing point w* of the distributions of the non equilibrium works. 

When the system is close to the critical point a definitive conclusion can not be stated. Data for small L suggest that 
when the correlation length approaches the system size L, the distribution p(Au) is very assymetric and there is no 
connection between the equilibrium value (Ait) and the crossing point w* of the non-equilibrium work distributions. A 
bigger computational effort (not affordable at present) would be needed in order to compute the accurate histograms 
for L > 24. The presented results may provide some clues for future investigations in order to extend work fluctuation 
theorems. 
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FIG. 6: (Color on line) Example of histograms corresponding to Pf(«o->i), Pr( — «i->o) an d P (Au — A(Hm)). Data corresponds 
to Ho — —W.0,Hi — —0.4, a = 3 and L = 12. Histograms are computed by cumulating 10 5 realizations of disorder. The 
arrows pointing down (up) indicate the peak (average value) of the distribution P(Au — A(HM)). 
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